Skip to content

Optimize calculation of 2-sample mean in MF_steps_ahead#83

Open
schristophe wants to merge 1 commit into
benfulcher:OperationChangesfrom
schristophe:main
Open

Optimize calculation of 2-sample mean in MF_steps_ahead#83
schristophe wants to merge 1 commit into
benfulcher:OperationChangesfrom
schristophe:main

Conversation

@schristophe

@schristophe schristophe commented Jun 18, 2026

Copy link
Copy Markdown

This PR optimises the calculation of the 2-sample mean in the MF_steps_ahead operation.

% A sliding mean of length 2
sm2p = zeros(N-i-1,1);
for j = 1:N-i-1
    seeds = yy(j:j+1);
    for k = 1:i % average with itself this many times
        p = mean(seeds);
        seeds = [seeds(2),p];
    end
    sm2p(j) = p;
end

Mathematically reframed, the code above iteratively calculates the terms of the following sequence (p subscript is for predict):

Let $y_j$ and $y_{j+1}$ be 2 consecutive samples of $y$.

$$y_\mathrm{p}(0) = y_j, y_\mathrm{p}(1) = y_{j+1},$$

$$\forall n\geq 2,\ y_\mathrm{p} (n) = \frac{y_\mathrm{p} (n-1) + y_\mathrm{p} (n-2)}{2},$$

which gives the value of the 2-sample mean predictor at $j+n$. This is a linear recurrence relation of order 2. It admits a closed-form expression that you can work out (see e.g. Wikipedia article).

$$\forall i\geq 1,\ y_p (i) = \frac{1}{3}\left(1+\frac{(-1)^{i+1}}{2^i}\right)y_j +\frac{2}{3}\left(1+\frac{(-1)^{i}}{2^{i+1}}\right)y_{j+1}$$

Note that I did a change of variable so that $i$ is the prediction length. Under this form, the "weights" used in the 2-sample mean can be computed directly regardless of the prediction length. This also implies that it is no longer needed to iterate "manually" over the samples of the time series using a for loop. The code above can then be rewritten as a simple matrix product:

% A sliding mean of length 2
weights = [1 + (-1)^(i+1)/2^i; 2 + (-1)^i/2^i] ./ 3;
sm2p = [yy(1:N-i-1), yy(2:N-i)] * weights;

This piece of code runs between tens to thousands of times faster on my laptop, depending on the length of the time series and the prediction length.
pr-hctsa-speedup

@schristophe
Optimize calculation of 2-sample mean in MF_steps_ahead
db0e417 
Previous implementation uses a 2nd order linear recurrence relation
to compute the 2-sample mean predictor at a given time index. It is
replaced by the closed-form solution that is faster to compute.
@benfulcher benfulcher changed the base branch from main to OperationChanges June 23, 2026 02:19
@benfulcher benfulcher changed the base branch from OperationChanges to main June 23, 2026 02:20
@benfulcher benfulcher changed the base branch from main to OperationChanges June 23, 2026 02:22
@benfulcher benfulcher changed the base branch from OperationChanges to main June 23, 2026 03:27
@benfulcher benfulcher changed the base branch from main to OperationChanges June 23, 2026 03:27
@benfulcher

benfulcher commented Jun 23, 2026

Copy link
Copy Markdown
Owner

Ok, thanks for the very clear PR! I assume you've validated this change in yielding equivalent performance? I don't see any evidence of this.
I think we indeed take the efficiency of this over the (potentially?) clearer for loop of what is actually being done. Have put it in the OperationChanges branch to minimize frequency of changes to main.

@schristophe

schristophe commented Jul 3, 2026 •
edited
Loading

Copy link
Copy Markdown
Author

I am not entirely sure what you mean about yielding equivalent performance. Perhaps I can add the results of some tests I did with the Bonn EGG dataset.

There are 4 master operations that are based on MF_steps_ahead in hctsa. All of them use a prediction length of 6. I computed them over the Bonn EGG dataset using the old code, then did the same with the "optimized" version. I used the timings in TS_CalcTime to do the comparison.

The first run of the MF_steps_ahead function (i.e. for the first time series and the first master operation) is notably slower in the optimized case (2.79 s) than in the old one (1.84 s). I only have basic knowledge of MATLAB, especially of how it works under the hood, so I am not sure what is causing this slowdown.

This is counterbalanced by all the other calls that are slightly faster on average. If I compare the timing for each time series this is what I get ($\Delta t = t_\mathrm{opt} - t_\mathrm{old}$):

$\Delta t$ MF_steps_ahead_ar_2_6 MF_steps_ahead_ar_best_6 MF_steps_ahead_ss_best_6 MF_steps_ahead_arma_3_1_6
mean -0.0370 -0.0353 -0.0398 -0.0369
std 0.0484 0.0185 0.0150 0.0184
min -0.0997 -0.1048 -0.0863 -0.1002
25% -0.0502 -0.0443 -0.0478 -0.0473
50% -0.0391 -0.0339 -0.0390 -0.0359
75% -0.0267 -0.0241 -0.0318 -0.0245
max 0.9561 0.0213 0.1295 0.0265

NB: One call to MF_steps_ahead() takes 0.1-0.3 s (roughly). Most of this time is spent fitting the model and calling predict()

For the whole dataset (N=500, 4 master operations), the speedup amounts to 74 s on a total computation time of ~6 min (18-19 s per master operation). I noticed that subsequent runs have a smaller speedup (0.025 s per time series, 10-12 s per master operation, 44 s for the whole dataset) because $t_\mathrm{old}$ is slightly lower on average. I do not know why.

Parallel computation is one scenario I have not tested.

I agree that the for loop makes it clearer to see how the "weights" are chosen. On the other hand, the matrix multiplication approach is more explicit about the values of the weights themselves and how they evolve with the prediction length. It is quite easy to see that the first weight oscillates around and converges to 1/3 and that the second has the same behavior but tends to 2/3.

I hope this helps. Let me know if there is anything else I can do.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

2 participants

@schristophe @benfulcher